USGS has libraries to retreive data from their system.

install libraries - do only one time when first runnign script

# https://github.com/USGS-R/dataRetrieval

# general R packages to install if you dont have them
# remove the "#" if you need to run the installations

# install.packages("devtools") # essential in installing other thigs
# install.packages("tidyverse") # installs a lot of things and ggplot
# install.packages("scales") # allows great scale formatting on ggplot
# install.packages("lubridate") # makes working with dates easier
# install.packages("janitor") # clean names of columns and other things
# install.packages("readxl") # read in excel files
# install.packages("plotly") # interactive plots
# install.packages(skimr)
# devtools::install_github("thomasp85/patchwork") # multiple plots

# specific to this module
# install.packages("dataRetrieval") # USGS Data Retreiveal Method

Libraires

# load these every time  
library(tidyverse)
library(scales)
library(lubridate)
library(janitor)
library(readxl)
library(skimr)
library(plotly)
library(patchwork)

# load the USGS package
library(dataRetrieval)

current data https://waterdata.usgs.gov/ny/nwis/uv/?site_no=01435000&PARAmeter_cd=00065,00060,63160

all data https://waterdata.usgs.gov/nwis/inventory/?site_no=01435000&agency_cd=USGS

Read in file for …


# The site is Neversink and the site number is 01435000
# USGS 01435000 NEVERSINK RIVER NEAR CLARYVILLE NY

site_no <-          "01435000"
par_code <-         "00060"
start.date <-      "1800-01-01"
end.date <-        "2050-01-01"

neversink_daily.df <- readNWISdv(siteNumbers = site_no,
                        parameterCd = par_code,
                        startDate = start.date,
                        endDate = end.date)
  
# rename the columns
neversink_daily.df <- renameNWISColumns(neversink_daily.df)

we need to make a date column to plot with

neversink_daily.df <- neversink_daily.df %>%
  select(agency_cd, 
         site_no, 
         date = Date, 
         discharge_cfs = Flow) 

we need to make a date column to plot with

neversink_daily.df <- neversink_daily.df %>%
  mutate(month = month(date))

Plot all data

q.plot <- neversink_daily.df %>% 
  ggplot(aes(date, discharge_cfs)) +
  geom_point(size=0.1) +
  geom_line()
q.plot

Plot February data

feb.plot <- neversink_daily.df %>% 
  filter(month == 2) %>%
  ggplot(aes(date, discharge_cfs)) +
  geom_point(size = 0.1) +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
feb.plot

NA

We can make it interactive

ggplotly(feb.plot)

August data

aug.plot <- neversink_daily.df %>%
  filter(month == 8) %>%
  ggplot(aes(date, discharge_cfs)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
aug.plot

Compare the plots

feb_aug.plot <- feb.plot +
                aug.plot +
                plot_layout(ncol =1)
feb_aug.plot

add trend lines as straight lines

feb_aug.plot <- feb.plot + geom_smooth(method="lm") +
                aug.plot + geom_smooth(method="lm") +
                plot_layout(ncol =1)
feb_aug.plot

How to see regression statistics

feb.model <- neversink_daily.df %>% 
  filter(month == 2) %>%
  lm(discharge_cfs ~ date, data=.)
summary(feb.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
-129.4  -86.5  -55.4    0.3 5925.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.585e+02  5.282e+00  30.007   <2e-16 ***
date        1.393e-03  5.713e-04   2.439   0.0148 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 233.5 on 2258 degrees of freedom
Multiple R-squared:  0.002627,  Adjusted R-squared:  0.002186 
F-statistic: 5.948 on 1 and 2258 DF,  p-value: 0.01481
aug.model <- neversink_daily.df %>% 
  filter(month == 8) %>%
  lm(discharge_cfs ~ date, data=.)
summary(aug.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
 -91.5  -61.5  -36.9   -2.6 7104.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.495e+01  4.794e+00  17.720  < 2e-16 ***
date        1.987e-03  5.152e-04   3.856 0.000118 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 220.2 on 2477 degrees of freedom
Multiple R-squared:  0.005966,  Adjusted R-squared:  0.005564 
F-statistic: 14.87 on 1 and 2477 DF,  p-value: 0.0001184
LS0tCnRpdGxlOiAiRWRkaWUgbW9kdWxlIHN0cmVhbSBkaXNjaGFyZ2UgZnJvbSBoYW5kb3V0IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKZWRpdG9yX29wdGlvbnM6IAogIGNodW5rX291dHB1dF90eXBlOiBpbmxpbmUKLS0tCgojIFVTR1MgaGFzIGxpYnJhcmllcyB0byByZXRyZWl2ZSBkYXRhIGZyb20gdGhlaXIgc3lzdGVtLiAgICAKCiMgaW5zdGFsbCBsaWJyYXJpZXMgLSBkbyBvbmx5IG9uZSB0aW1lIHdoZW4gZmlyc3QgcnVubmlnbiBzY3JpcHQgICAgCmBgYHtyIGluc3RhbGwgcGFja2FnZXMsIG1lc3NhZ2U9VFJVRSwgd2FybmluZz1UUlVFfQojIGh0dHBzOi8vZ2l0aHViLmNvbS9VU0dTLVIvZGF0YVJldHJpZXZhbAoKIyBnZW5lcmFsIFIgcGFja2FnZXMgdG8gaW5zdGFsbCBpZiB5b3UgZG9udCBoYXZlIHRoZW0KIyByZW1vdmUgdGhlICIjIiBpZiB5b3UgbmVlZCB0byBydW4gdGhlIGluc3RhbGxhdGlvbnMKCiMgaW5zdGFsbC5wYWNrYWdlcygiZGV2dG9vbHMiKSAjIGVzc2VudGlhbCBpbiBpbnN0YWxsaW5nIG90aGVyIHRoaWdzCiMgaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIikgIyBpbnN0YWxscyBhIGxvdCBvZiB0aGluZ3MgYW5kIGdncGxvdAojIGluc3RhbGwucGFja2FnZXMoInNjYWxlcyIpICMgYWxsb3dzIGdyZWF0IHNjYWxlIGZvcm1hdHRpbmcgb24gZ2dwbG90CiMgaW5zdGFsbC5wYWNrYWdlcygibHVicmlkYXRlIikgIyBtYWtlcyB3b3JraW5nIHdpdGggZGF0ZXMgZWFzaWVyCiMgaW5zdGFsbC5wYWNrYWdlcygiamFuaXRvciIpICMgY2xlYW4gbmFtZXMgb2YgY29sdW1ucyBhbmQgb3RoZXIgdGhpbmdzCiMgaW5zdGFsbC5wYWNrYWdlcygicmVhZHhsIikgIyByZWFkIGluIGV4Y2VsIGZpbGVzCiMgaW5zdGFsbC5wYWNrYWdlcygicGxvdGx5IikgIyBpbnRlcmFjdGl2ZSBwbG90cwojIGluc3RhbGwucGFja2FnZXMoc2tpbXIpCiMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJ0aG9tYXNwODUvcGF0Y2h3b3JrIikgIyBtdWx0aXBsZSBwbG90cwoKIyBzcGVjaWZpYyB0byB0aGlzIG1vZHVsZQojIGluc3RhbGwucGFja2FnZXMoImRhdGFSZXRyaWV2YWwiKSAjIFVTR1MgRGF0YSBSZXRyZWl2ZWFsIE1ldGhvZApgYGAKCgpMaWJyYWlyZXMKYGBge3IgbGFvZCBsaWJyYXJpZXMsIG1lc3NhZ2U9VFJVRSwgd2FybmluZz1UUlVFfQojIGxvYWQgdGhlc2UgZXZlcnkgdGltZSAgCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHNjYWxlcykKbGlicmFyeShsdWJyaWRhdGUpCmxpYnJhcnkoamFuaXRvcikKbGlicmFyeShyZWFkeGwpCmxpYnJhcnkoc2tpbXIpCmxpYnJhcnkocGxvdGx5KQpsaWJyYXJ5KHBhdGNod29yaykKCiMgbG9hZCB0aGUgVVNHUyBwYWNrYWdlCmxpYnJhcnkoZGF0YVJldHJpZXZhbCkKYGBgCgpjdXJyZW50IGRhdGEKaHR0cHM6Ly93YXRlcmRhdGEudXNncy5nb3Yvbnkvbndpcy91di8/c2l0ZV9ubz0wMTQzNTAwMCZQQVJBbWV0ZXJfY2Q9MDAwNjUsMDAwNjAsNjMxNjAKCmFsbCBkYXRhCmh0dHBzOi8vd2F0ZXJkYXRhLnVzZ3MuZ292L253aXMvaW52ZW50b3J5Lz9zaXRlX25vPTAxNDM1MDAwJmFnZW5jeV9jZD1VU0dTCgoKUmVhZCBpbiBmaWxlIGZvciAuLi4KYGBge3IgcmVhZCBpbiBmaWxlfQoKIyBUaGUgc2l0ZSBpcyBOZXZlcnNpbmsgYW5kIHRoZSBzaXRlIG51bWJlciBpcyAwMTQzNTAwMAojIFVTR1MgMDE0MzUwMDAgTkVWRVJTSU5LIFJJVkVSIE5FQVIgQ0xBUllWSUxMRSBOWQoKc2l0ZV9ubyA8LSAgICAgICAgICAiMDE0MzUwMDAiCnBhcl9jb2RlIDwtICAgICAgICAgIjAwMDYwIgpzdGFydC5kYXRlIDwtICAgICAgIjE4MDAtMDEtMDEiCmVuZC5kYXRlIDwtICAgICAgICAiMjA1MC0wMS0wMSIKCm5ldmVyc2lua19kYWlseS5kZiA8LSByZWFkTldJU2R2KHNpdGVOdW1iZXJzID0gc2l0ZV9ubywKICAgICAgICAgICAgICAgICAgICAgICAgcGFyYW1ldGVyQ2QgPSBwYXJfY29kZSwKICAgICAgICAgICAgICAgICAgICAgICAgc3RhcnREYXRlID0gc3RhcnQuZGF0ZSwKICAgICAgICAgICAgICAgICAgICAgICAgZW5kRGF0ZSA9IGVuZC5kYXRlKQogIAojIHJlbmFtZSB0aGUgY29sdW1ucwpuZXZlcnNpbmtfZGFpbHkuZGYgPC0gcmVuYW1lTldJU0NvbHVtbnMobmV2ZXJzaW5rX2RhaWx5LmRmKQpgYGAKCiMgd2UgbmVlZCB0byBtYWtlIGEgZGF0ZSBjb2x1bW4gdG8gcGxvdCB3aXRoCmBgYHtyfQpuZXZlcnNpbmtfZGFpbHkuZGYgPC0gbmV2ZXJzaW5rX2RhaWx5LmRmICU+JQogIHNlbGVjdChhZ2VuY3lfY2QsIAogICAgICAgICBzaXRlX25vLCAKICAgICAgICAgZGF0ZSA9IERhdGUsIAogICAgICAgICBkaXNjaGFyZ2VfY2ZzID0gRmxvdykgCmBgYAoKIyB3ZSBuZWVkIHRvIG1ha2UgYSBkYXRlIGNvbHVtbiB0byBwbG90IHdpdGgKYGBge3J9Cm5ldmVyc2lua19kYWlseS5kZiA8LSBuZXZlcnNpbmtfZGFpbHkuZGYgJT4lCiAgbXV0YXRlKG1vbnRoID0gbW9udGgoZGF0ZSkpCmBgYAoKCgojIFBsb3QgYWxsIGRhdGEKYGBge3J9CnEucGxvdCA8LSBuZXZlcnNpbmtfZGFpbHkuZGYgJT4lIAogIGdncGxvdChhZXMoZGF0ZSwgZGlzY2hhcmdlX2NmcykpICsKICBnZW9tX3BvaW50KHNpemU9MC4xKSArCiAgZ2VvbV9saW5lKCkKcS5wbG90CmBgYAoKIyBQbG90IEZlYnJ1YXJ5IGRhdGEgCmBgYHtyfQpmZWIucGxvdCA8LSBuZXZlcnNpbmtfZGFpbHkuZGYgJT4lIAogIGZpbHRlcihtb250aCA9PSAyKSAlPiUKICBnZ3Bsb3QoYWVzKGRhdGUsIGRpc2NoYXJnZV9jZnMpKSArCiAgZ2VvbV9wb2ludChzaXplID0gMC4xKSArCiAgZ2VvbV9saW5lKCkgKwogIGxhYnMoeT0iTWVhbiBNb250aGx5IERpc2NoYXJnZSBtXjMvc2VjIiwgeD0iRGF0ZSIpCmZlYi5wbG90CiAgCmBgYAoKV2UgY2FuIG1ha2UgaXQgaW50ZXJhY3RpdmUKYGBge3J9CmdncGxvdGx5KGZlYi5wbG90KQpgYGAKCgpBdWd1c3QgZGF0YQpgYGB7cn0KYXVnLnBsb3QgPC0gbmV2ZXJzaW5rX2RhaWx5LmRmICU+JQogIGZpbHRlcihtb250aCA9PSA4KSAlPiUKICBnZ3Bsb3QoYWVzKGRhdGUsIGRpc2NoYXJnZV9jZnMpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2xpbmUoKSArCiAgbGFicyh5PSJNZWFuIE1vbnRobHkgRGlzY2hhcmdlIG1eMy9zZWMiLCB4PSJEYXRlIikKYXVnLnBsb3QKYGBgCgojIENvbXBhcmUgdGhlIHBsb3RzCmBgYHtyfQpmZWJfYXVnLnBsb3QgPC0gZmViLnBsb3QgKwogICAgICAgICAgICAgICAgYXVnLnBsb3QgKwogICAgICAgICAgICAgICAgcGxvdF9sYXlvdXQobmNvbCA9MSkKZmViX2F1Zy5wbG90CmBgYAoKIyBhZGQgdHJlbmQgbGluZXMgYXMgc3RyYWlnaHQgbGluZXMKYGBge3J9CmZlYl9hdWcucGxvdCA8LSBmZWIucGxvdCArIGdlb21fc21vb3RoKG1ldGhvZD0ibG0iKSArCiAgICAgICAgICAgICAgICBhdWcucGxvdCArIGdlb21fc21vb3RoKG1ldGhvZD0ibG0iKSArCiAgICAgICAgICAgICAgICBwbG90X2xheW91dChuY29sID0xKQpmZWJfYXVnLnBsb3QKYGBgCgojIEhvdyB0byBzZWUgcmVncmVzc2lvbiBzdGF0aXN0aWNzCmBgYHtyfQpmZWIubW9kZWwgPC0gbmV2ZXJzaW5rX2RhaWx5LmRmICU+JSAKICBmaWx0ZXIobW9udGggPT0gMikgJT4lCiAgbG0oZGlzY2hhcmdlX2NmcyB+IGRhdGUsIGRhdGE9LikKc3VtbWFyeShmZWIubW9kZWwpCgpgYGAKCgpgYGB7cn0KYXVnLm1vZGVsIDwtIG5ldmVyc2lua19kYWlseS5kZiAlPiUgCiAgZmlsdGVyKG1vbnRoID09IDgpICU+JQogIGxtKGRpc2NoYXJnZV9jZnMgfiBkYXRlLCBkYXRhPS4pCnN1bW1hcnkoYXVnLm1vZGVsKQoKYGBgCgoK